resultsPath=file.path(getwd(),"Results")
# Gather parameters from command line
#dir.create(file.path(resultsPath,"cache"), showWarnings=F, recursive=T)
nCores <- params$nCores #parallel::detectCores()
subsetGenes <- params$subsetGenes
subsetCells <- params$subsetCells
resolution <- as.numeric(params$resolution)
interactive <- params$interactive
perplexity <- params$perplexity
perplexity <- ifelse(exists("perplexity"), perplexity, 30)
root <- getwd()
# Have to setwd via knitr
# knitr::opts_knit$set(root.dir=resultsPath, child.path = resultsPath)
knitr::opts_chunk$set(echo=T, error=T, root.dir = resultsPath
# cache=T, cache.lazy=T
)
# Utilize parallel processing later on
cat(paste("**** Utilized Cores **** =", nCores)) ## **** Utilized Cores **** = 2
params## $subsetGenes
## [1] "protein_coding"
##
## $subsetCells
## [1] 500
##
## $resolution
## [1] 0.6
##
## $resultsPath
## [1] "./"
##
## $nCores
## [1] 2
##
## $perplexity
## [1] 30
** ./ **
library(Seurat)
library(dplyr)
library(gridExtra)
library(knitr)
library(plotly)
library(ggplot2)
library(viridis)
library(reshape2)
library(shiny)
library(ggrepel)
library(DT)
library(ComplexHeatmap); #BiocManager::install("ComplexHeatmap")
# install.packages('devtools')
# devtools::install_github('talgalili/heatmaply')
## Install Bioconductor
# if (!requireNamespace("BiocManager"))
# install.packages("BiocManager")
# BiocManager::install(c("biomaRt"))
library(biomaRt)
# BiocManager::install(c("DESeq2"))
library(DESeq2)
# library(snow); #BiocManager::install("Rmpi") #NOTE: different lib name than install name (snow vs Rmpi)
library(enrichR) #BiocManager::install("enrichR")
# Useful Seurat functions
## Seurat::FindGeneTerms() # Enrichr API
## Seurat::MultiModal_CCA() # Integrates data from disparate datasets (CIA version too)
sessionInfo()## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.14.2
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 grid stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] enrichR_1.0 DESeq2_1.22.2
## [3] SummarizedExperiment_1.12.0 DelayedArray_0.8.0
## [5] BiocParallel_1.16.5 matrixStats_0.54.0
## [7] Biobase_2.42.0 GenomicRanges_1.34.0
## [9] GenomeInfoDb_1.18.1 IRanges_2.16.0
## [11] S4Vectors_0.20.1 BiocGenerics_0.28.0
## [13] biomaRt_2.38.0 ComplexHeatmap_1.20.0
## [15] DT_0.5.1 ggrepel_0.8.0
## [17] shiny_1.2.0 reshape2_1.4.3
## [19] viridis_0.5.1 viridisLite_0.3.0
## [21] plotly_4.8.0 knitr_1.21
## [23] gridExtra_2.3 dplyr_0.7.8
## [25] Seurat_2.3.4 Matrix_1.2-15
## [27] cowplot_0.9.4 ggplot2_3.1.0
##
## loaded via a namespace (and not attached):
## [1] snow_0.4-3 backports_1.1.3 circlize_0.4.5
## [4] Hmisc_4.1-1 plyr_1.8.4 igraph_1.2.2
## [7] lazyeval_0.2.1 splines_3.5.1 digest_0.6.18
## [10] foreach_1.4.4 htmltools_0.3.6 lars_1.2
## [13] gdata_2.18.0 memoise_1.1.0 magrittr_1.5
## [16] checkmate_1.9.1 cluster_2.0.7-1 mixtools_1.1.0
## [19] ROCR_1.0-7 annotate_1.60.0 R.utils_2.7.0
## [22] prettyunits_1.0.2 colorspace_1.4-0 blob_1.1.1
## [25] xfun_0.4 crayon_1.3.4 RCurl_1.95-4.11
## [28] jsonlite_1.6 genefilter_1.64.0 bindr_0.1.1
## [31] survival_2.43-3 zoo_1.8-4 iterators_1.0.10
## [34] ape_5.2 glue_1.3.0 gtable_0.2.0
## [37] zlibbioc_1.28.0 XVector_0.22.0 GetoptLong_0.1.7
## [40] kernlab_0.9-27 shape_1.4.4 prabclus_2.2-7
## [43] DEoptimR_1.0-8 scales_1.0.0 mvtnorm_1.0-8
## [46] DBI_1.0.0 bibtex_0.4.2 Rcpp_1.0.0
## [49] metap_1.0 dtw_1.20-1 progress_1.2.0
## [52] xtable_1.8-3 htmlTable_1.13.1 reticulate_1.10
## [55] foreign_0.8-71 bit_1.1-14 proxy_0.4-22
## [58] mclust_5.4.2 SDMTools_1.1-221 Formula_1.2-3
## [61] tsne_0.1-3 htmlwidgets_1.3 httr_1.4.0
## [64] gplots_3.0.1 RColorBrewer_1.1-2 fpc_2.1-11.1
## [67] acepack_1.4.1 modeltools_0.2-22 ica_1.0-2
## [70] pkgconfig_2.0.2 XML_3.98-1.16 R.methodsS3_1.7.1
## [73] flexmix_2.3-14 nnet_7.3-12 locfit_1.5-9.1
## [76] tidyselect_0.2.5 rlang_0.3.1 later_0.7.5
## [79] AnnotationDbi_1.44.0 munsell_0.5.0 tools_3.5.1
## [82] RSQLite_2.1.1 ggridges_0.5.1 evaluate_0.12
## [85] stringr_1.3.1 yaml_2.2.0 npsurv_0.4-0
## [88] bit64_0.9-7 fitdistrplus_1.0-11 robustbase_0.93-3
## [91] caTools_1.17.1.1 purrr_0.2.5 RANN_2.6.1
## [94] bindrcpp_0.2.2 pbapply_1.3-4 nlme_3.1-137
## [97] mime_0.6 R.oo_1.22.0 hdf5r_1.0.1
## [100] compiler_3.5.1 rstudioapi_0.9.0 png_0.1-7
## [103] lsei_1.2-0 geneplotter_1.60.0 tibble_2.0.1
## [106] stringi_1.2.4 lattice_0.20-38 trimcluster_0.1-2.1
## [109] pillar_1.3.1 Rdpack_0.10-1 lmtest_0.9-36
## [112] GlobalOptions_0.1.0 data.table_1.12.0 bitops_1.0-6
## [115] irlba_2.3.2 gbRd_0.4-11 httpuv_1.4.5.1
## [118] R6_2.3.0 latticeExtra_0.6-28 promises_1.0.1
## [121] KernSmooth_2.23-15 codetools_0.2-16 MASS_7.3-51.1
## [124] gtools_3.8.1 assertthat_0.2.0 rjson_0.2.20
## [127] withr_2.1.2 GenomeInfoDbData_1.2.0 hms_0.4.2
## [130] diptest_0.75-7 doSNOW_1.0.16 rpart_4.1-13
## [133] tidyr_0.8.2 class_7.3-15 rmarkdown_1.11
## [136] segmented_0.5-3.0 Rtsne_0.15 base64enc_0.1-3
print(paste("Seurat ", packageVersion("Seurat")))## [1] "Seurat 2.3.4"
createDT <- function(DF, caption="", scrollY=500){
data <- DT::datatable(DF, caption=caption,
extensions = 'Buttons',
options = list( dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
scrollY = scrollY, scrollX=T, scrollCollapse = T, paging = F,
columnDefs = list(list(className = 'dt-center', targets = "_all"))
)
)
return(data)
}Rstudio has a default memory limit of only 1GB. To override this, detect the true memory available and set a new limit.
library(ulimit) # devtools::install_github("krlmlr/ulimit")
library(benchmarkme)
RAM <- print(benchmarkme::get_ram())
## Convert GB to Mib
RAM_Mib <- strsplit(RAM, " ")[[1]][1] %>% as.numeric() * 953.67431640625
cat(paste("Available RAM:",RAM))
## Set new memory limit
ulimit::memory_limit(RAM_Mib) ## ! IMPORTANT! Must not setwd to local path when launching on cluster
# setwd("~/Desktop/PD_scRNAseq/")
dir.create(file.path(root,"Data"), showWarnings=F)
load(file.path(root,"Data/seurat_object_add_HTO_ids.Rdata"))
DAT <- seurat.obj
rm(seurat.obj)DAT## An object of class seurat in project RAJ_13357
## 24914 genes across 22113 samples.
metadata <- read.table(file.path(root,"Data/meta.data4.tsv"))
createDT( metadata, caption = "Metadata") ## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html
# Make AgeGroups
makeAgeGroups <- function(){
dim(metadata)
getMaxRound <- function(vals=metadata$Age, unit=10)unit*ceiling((max(vals)/unit))
getMinRound <- function(vals=metadata$Age, unit=10)unit*floor((min(vals)/unit))
ageBreaks = c(seq(getMinRound(), getMaxRound(), by = 10), getMaxRound()+10)
AgeGroupsUniq <- c()
for (i in 1:(length(ageBreaks)-1)){
AgeGroupsUniq <- append(AgeGroupsUniq, paste(ageBreaks[i],ageBreaks[i+1], sep="-"))
}
data.table::setDT(metadata,keep.rownames = T,check.names = F)[, AgeGroups := cut(Age,
breaks = ageBreaks,
right = F,
labels = AgeGroupsUniq,
nclude.lowest=T)]
metadata <- data.frame(metadata)
unique(metadata$AgeGroups)
head(metadata)
dim(metadata)
return(metadata)
}
# metadata <- makeAgeGroups()
DAT <- AddMetaData(object = DAT, metadata = metadata)
# Get rid of any NAs (cells that don't match up with the metadata)
if(subsetCells==F){
DAT <- FilterCells(object = DAT, subset.names = "nGene", low.thresholds = 0)
} else {DAT <- FilterCells(object = DAT, subset.names = "nGene", low.thresholds = 0,
# Subset for testing
cells.use = DAT@cell.names[0:subsetCells]
)
} Include only subsets of genes by type. Biotypes from: https://useast.ensembl.org/info/genome/genebuild/biotypes.html
subsetBiotypes <- function(DAT, subsetGenes){
if( subsetGenes!=F ){
cat(paste("Subsetting genes:",subsetGenes, "\n"))
# If the gene_biotypes file exists, import csv. Otherwise, get from biomaRt
if(file_test("-f", file.path(root,"Data/gene_biotypes.csv"))){
biotypes <- read.csv(file.path(root,"Data/gene_biotypes.csv"))
}
else {
ensembl <- useMart(biomart="ENSEMBL_MART_ENSEMBL", host="grch37.ensembl.org",
dataset="hsapiens_gene_ensembl")
ensembl <- useDataset(mart = ensembl, dataset = "hsapiens_gene_ensembl")
listFilters(ensembl)
listAttributes(ensembl)
biotypes <- getBM(attributes=c("hgnc_symbol", "gene_biotype"), filters="hgnc_symbol",
values=row.names(DAT@data), mart=ensembl)
write.csv(biotypes, file.path(root,"Data/gene_biotypes.csv"), quote=F, row.names=F)
}
# Subset data by creating new Seurat object (annoying but necessary)
geneSubset <- biotypes[biotypes$gene_biotype==subsetGenes,"hgnc_symbol"]
cat(paste(dim(DAT@raw.data[geneSubset, ])[1],"/", dim(DAT@raw.data)[1],
"genes are", subsetGenes))
# Add back into DAT
subset.matrix <- DAT@raw.data[geneSubset, ] # Pull the raw expression matrix from the original Seurat object containing only the genes of interest
DAT_sub <- CreateSeuratObject(subset.matrix) # Create a new Seurat object with just the genes of interest
orig.ident <- row.names(DAT@meta.data) # Pull the identities from the original Seurat object as a data.frame
DAT_sub <- AddMetaData(object = DAT_sub, metadata = DAT@meta.data) # Add the idents to the meta.data slot
DAT_sub <- SetAllIdent(object = DAT_sub, id = "ident") # Assign identities for the new Seurat object
DAT <- DAT_sub
rm(list = c("DAT_sub","geneSubset", "subset.matrix", "orig.ident"))
}
}
subsetBiotypes(DAT, subsetGenes)## Subsetting genes: protein_coding
## 14827 / 24914 genes are protein_coding
Filter by cells, normalize , filter by gene variability.
cat("Total Cells:", length(DAT@cell.names), "\n")## Total Cells: 495
DAT <- FilterCells(object = DAT, subset.names = c("nGene", "percent.mito"),
low.thresholds = c(200, -Inf), high.thresholds = c(2500, 0.05))
cat("Filtered Cells:", length(DAT@cell.names))## Filtered Cells: 495
DAT <- NormalizeData(object = DAT, normalization.method = "LogNormalize",
scale.factor = 10000)** Important!**: * In ScaleData… + Specify do.par = F (unless you have parallel processing set up properly, this will cause your script to crash) + Specify num.cores = nCores (to use all available cores, determined by parallel::detectCores())
Regress out: number of unique transcripts (nUMI), % mitochondrial transcripts (percent.mito)
# Store the top most variable genes in @var.genes
DAT <- FindVariableGenes(object = DAT, mean.function = ExpMean, dispersion.function = LogVMR,
x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)cat("Total Genes:", length(row.names(DAT@raw.data)), "\n")## Total Genes: 24914
cat("Highly Variable Genes:", length(DAT@var.genes))## Highly Variable Genes: 9672
# IMPORTANT!: Must set do.par=T and num.cors = n for large datasets being processed on computing clusters
# IMPORTANT!: Use only the var.genes identified by 'FindVariableGenes' as the 'gene.use' arg in 'ScaleData'
## This will greatly reduced the computational load.
# Ensure CD14 and CD16 are included
appendedGenes <- c(DAT@var.genes, "CD14", "FCGR3A")
DAT <- ScaleData(object = DAT, genes.use = appendedGenes , vars.to.regress = c("nUMI", "percent.mito"),
do.par = T, num.cores = nCores)## Regressing out: nUMI, percent.mito
##
## Time Elapsed: 5.45868277549744 secs
## Scaling data matrix
DAT## An object of class seurat in project RAJ_13357
## 24914 genes across 495 samples.
vp <- VlnPlot(object = DAT, features.plot = c("nGene", "nUMI", "percent.mito"),nCol = 3, do.return = T) %>% + ggplot2::aes(alpha=0.5)
vp# par(mfrow = c(1, 2))
# do.hover <- ifelse(interactive==T, T, F)
gp1 <- GenePlot(object = DAT, gene1 = "nUMI", gene2 = "percent.mito", pch.use=20) #do.hover=do.hover, data.hover = "mut")gp1## NULL
# do.hover <-ifelse(interactive==T, T, F)
gp2 <- GenePlot(object = DAT, gene1 = "nUMI", gene2 = "nGene", pch.use=20) #do.hover=do.hover, data.hover = "mut")gp2## NULL
ProjectPCA scores each gene in the dataset (including genes not included in the PCA) based on their correlation with the calculated components. Though we don’t use this further here, it can be used to identify markers that are strongly correlated with cellular heterogeneity, but may not have passed through variable gene selection. The results of the projected PCA can be explored by setting use.full=T in the functions above
# Run PCA with only the top most variables genes
DAT <- RunPCA(object = DAT, pc.genes = DAT@var.genes, do.print=F, verbose=F) #, pcs.print = 1:5, genes.print = 5
# Store in Seurat object so you don't have to recalculate it for the tSNE/UMAP steps
DAT <- ProjectPCA(object = DAT, do.print=F) VizPCA(object = DAT, pcs.use = 1:2)# 'PCHeatmap' is a wrapper for heatmap.2
PCHeatmap(object = DAT, pc.use = 1:12,do.balanced=T, label.columns=F, use.full=F) # cells.use = 500, Determine statistically significant PCs for further analysis. NOTE: This process can take a long time for big datasets, comment out for expediency. More approximate techniques such as those implemented in PCElbowPlot() can be used to reduce computation time
#DAT <- JackStraw(object = DAT, num.replicate = 100, display.progress = FALSE)
PCElbowPlot(object = DAT)We first construct a KNN graph based on the euclidean distance in PCA space, and refine the edge weights between any two cells based on the shared overlap in their local neighborhoods (Jaccard similarity). To cluster the cells, we apply modularity optimization techniques such as the Louvain algorithm (default) or SLM [SLM, Blondel et al., Journal of Statistical Mechanics], to iteratively group cells together, with the goal of optimizing the standard modularity function.
On Resolution
The FindClusters function implements the procedure, and contains a resolution parameter that sets the ‘granularity’ of the downstream clustering, with increased values leading to a greater number of clusters. We find that setting this parameter between 0.6-1.2 typically returns good results for single cell datasets of around 3K cells. Optimal resolution often increases for larger datasets. The clusters are saved in the object@ident slot.
The clustering approach in FindClusters was heavily inspired by recent manuscripts which applied graph-based clustering approaches to scRNA-seq data SNN-Cliq, Xu and Su, Bioinformatics, 2015 and CyTOF data PhenoGraph, Levine et al., Cell, 2015.
DAT <- RunTSNE(object=DAT, reduction.use = "pca", dims.use = 1:10, do.fast = TRUE,
perplexity = perplexity, tsne.method = "Rtsne", num_threads=nCores, verbose=F) # FItSNE
# TRY DIFFERENT RESOLUTIONS
DAT <- StashIdent(object = DAT, save.name = "pre_clustering")
# DAT <- SetAllIdent(object = DAT, id = "pre_clustering")
DAT <- FindClusters(object = DAT, reduction.type = "pca", dims.use = 1:10, algorithm = 1,
resolution = resolution, print.output = T, save.SNN = T,
n.start = 10, nn.eps = 0.5, plot.SNN = T, force.recalc=T) ## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 495
## Number of edges: 37091
##
## Running Louvain algorithm...
## Random start: 1
## Iteration: 1
## Modularity: 0.5038
## Iteration: 2
## Modularity: 0.5038
##
## Random start: 2
## Iteration: 1
## Modularity: 0.5038
## Iteration: 2
## Modularity: 0.5038
##
## Random start: 3
## Iteration: 1
## Modularity: 0.5038
## Iteration: 2
## Modularity: 0.5038
##
## Random start: 4
## Iteration: 1
## Modularity: 0.5038
## Iteration: 2
## Modularity: 0.5038
##
## Random start: 5
## Iteration: 1
## Modularity: 0.5038
## Iteration: 2
## Modularity: 0.5038
##
## Random start: 6
## Iteration: 1
## Modularity: 0.5038
## Iteration: 2
## Modularity: 0.5038
##
## Random start: 7
## Iteration: 1
## Modularity: 0.5020
## Iteration: 2
## Modularity: 0.5039
## Iteration: 3
## Modularity: 0.5039
##
## Random start: 8
## Iteration: 1
## Modularity: 0.5038
## Iteration: 2
## Modularity: 0.5038
##
## Random start: 9
## Iteration: 1
## Modularity: 0.5039
## Iteration: 2
## Modularity: 0.5039
##
## Random start: 10
## Iteration: 1
## Modularity: 0.5038
## Iteration: 2
## Modularity: 0.5038
##
## Maximum modularity in 10 random starts: 0.5039
## Number of communities: 3
## Elapsed time: 0 seconds
PrintFindClustersParams(object = DAT) ## Parameters used in latest FindClusters calculation run on: 2019-01-22 20:19:59
## =============================================================================
## Resolution: 0.6
## -----------------------------------------------------------------------------
## Modularity Function Algorithm n.start n.iter
## 1 1 10 10
## -----------------------------------------------------------------------------
## Reduction used k.param prune.SNN
## pca 30 0.0667
## -----------------------------------------------------------------------------
## Dims used in calculation
## =============================================================================
## 1 2 3 4 5 6 7 8 9 10
DAT <- StashIdent(object = DAT, save.name = "post_clustering") # do.hover <-ifelse(interactive==T, T, F)
PCAPlot(object = DAT, dim.1 = 1, dim.2 = 2, group.by="post_clustering")#, do.hover=do.hover, data.hover="mut")Additional UMAP arguments detailed here: https://umap-learn.readthedocs.io/en/latest/api.html#module-umap.umap_
# cat(print(mem_used()))
DAT <- RunUMAP(object = DAT, reduction.use = "pca", dims.use = 1:10, verbose=TRUE, num_threads=nCores) # , num_threads=0
# cat(print(mem_used()))
# Plot results
DimPlot(object = DAT, reduction.use = 'umap')As input to the tSNE, we suggest using the same PCs as input to the clustering analysis, although computing the tSNE based on scaled gene expression is also supported using the genes.use argument.
** Important!**: Specify num_threads=0 in ‘RunTSNE’ to use all available cores.
“FItSNE”, a new fast implementation of t-SNE, is also available through RunTSNE. However FItSNE must first be setup on your computer.
labSize <- 12
customColors <- function(var="post_clustering", palette="Set1"){
add.alpha <- function(col, alpha=1){
if(missing(col))
stop("Please provide a vector of colours.")
apply(sapply(col, col2rgb)/255, 2,
function(x)
rgb(x[1], x[2], x[3], alpha=alpha))
}
cluster_colors <- RColorBrewer::brewer.pal( length(unique(DAT@meta.data[var])), palette)
cluster_colors_transparent <- add.alpha(cluster_colors, .5) %>% as.character()
return(cluster_colors_transparent)
}
# Try t-SNE at different perplexities
for (i in c(perplexity,5,20,30,100)){
cat('\n')
cat("### t-SNE: perplexity =",i,"\n")
DAT <- RunTSNE(object=DAT, reduction.use = "pca", dims.use = 1:10, do.fast = TRUE,
perplexity = i, tsne.method = "Rtsne", num_threads=nCores, verbose=F) # FItSNE
tsnePlot <- TSNEPlot(object = DAT, do.label=T, label.size = labSize, do.return=T) +
scale_color_brewer( palette = "Set1", aesthetics = aes(alpha=.5))
print(tsnePlot)
cat('\n')
} tSNE_metadata_plot <- function(var, labSize=12){
cat(paste("t-SNE Metadata plot for ", var))
# Metadata plot
p1 <- TSNEPlot(DAT, do.return = T, do.label = T, group.by = var,label.size = labSize,
plot.title=paste(var), vector.friendly=T) +
theme(legend.position = "top") +
scale_color_brewer( palette = "Dark2", aesthetics = aes(alpha=.5))
# t-SNE clusters plot
p2 <- TSNEPlot(DAT, do.return = T, do.label = T,label.size = labSize,
plot.title=paste("Unsupervised Clusters"), vector.friendly=T) +
theme(legend.position = "top") +
scale_color_brewer( palette = "Set1", aesthetics = aes(alpha=.5))
print(plot_grid(p1,p2))
}
# Iterate plots over metadata variables
metaVars <- c("dx","mut","Gender","Age", "ID")
for (var in metaVars){
cat('\n')
cat("### t-SNE Metadata plot for ",var, "\n")
tSNE_metadata_plot(var)
cat('\n')
} ##
## ### t-SNE Metadata plot for dx
## t-SNE Metadata plot for dx
##
##
## ### t-SNE Metadata plot for mut
## t-SNE Metadata plot for mut
##
##
## ### t-SNE Metadata plot for Gender
## t-SNE Metadata plot for Gender
##
##
## ### t-SNE Metadata plot for Age
## t-SNE Metadata plot for Age
##
##
## ### t-SNE Metadata plot for ID
## t-SNE Metadata plot for ID
Seurat has several tests for differential expression which can be set with the test.use parameter (see the DE vignette for details). For example, the ROC test returns the ‘classification power’ for any individual marker (ranging from 0 - random, to 1 - perfect).
Shown here: Biomarkers of each cluster vs. all other clusters.
# Limit to only top variable genes?:
# Set arg 'only.pos=F' to capture negative biomarkers
# DAT <- SetIdent(DAT, ident.use = "post_clustering")
DAT.markers <- FindAllMarkers(object = DAT, min.pct = 0.25, thresh.use = 0.25, only.pos = F, test.use = "wilcox")
DAT.markers <- DAT.markers %>% mutate(FC = 2^avg_logFC)
DAT.markers.sig <- DAT.markers %>% subset(p_val_adj<=0.05)
markers.summary <- DAT.markers.sig %>% group_by(cluster) %>% tally()
# markers.summary <- base::merge(DAT.markers.sig %>% group_by(cluster) %>% tally(),
# DAT.markers %>% group_by(cluster) %>% summarise(mean(avg_logFC)),
# by="cluster" )
createDT(markers.summary, caption = "Number of DEGs and Mean logFC per Cluster")createDT(DAT.markers, caption = paste("All Biomarkers: All Clusters"))topNum = 5
topBiomarkers <- DAT.markers %>% group_by(cluster) %>% top_n(topNum, avg_logFC)
createDT(DAT.markers, caption = paste("All Biomarkers: All Clusters"))getTopBiomarker <- function(DAT.markers, clusterID, topN=1){
df <-DAT.markers %>%
subset(p_val_adj<0.05 & cluster==as.character(clusterID)) %>%
arrange(desc(avg_logFC))
top_pct_markers <- df[1:topN,"gene"]
return(top_pct_markers)
}
# clust1_biomarkers <- getTopBiomarker(DAT.markers, clusterID=1, topN=2)
# clust2_biomarkers <- getTopBiomarker(DAT.markers, clusterID=2, topN=2)
### Plot biomarkers
plotBiomarkers <- function(DAT, biomarkers, cluster){
biomarkerPlots <- list()
for (marker in biomarkers){
p <- VlnPlot(object = DAT, features.plot = c(marker), y.log=T, return.plotlist=T)
biomarkerPlots[[marker]] <- p + ggplot2::aes(alpha=0.5) + xlab( "Cluster") + ylab( "Expression")
}
combinedPlot <- do.call(grid.arrange, c(biomarkerPlots, list(ncol=2, top=paste("Top DEG Biomarkers for Cluster",cluster))) )
# biomarkerPlots <- lapply(biomarkers, function(marker) {
# VlnPlot(object = DAT, features.plot = c(marker), y.log=T, return.plotlist=T) %>% + ggplot2::ggtitle(marker) %>% ggplotly()
# })
# return(subplot(biomarkerPlots) )
}
top1 <- DAT.markers %>% group_by(cluster) %>% top_n(1, avg_logFC)
nCols <- floor( sqrt(length(unique(top1$cluster))) )
figHeight <- nCols *7
# Plot top 2 biomarker genes for each
for (clust in unique(DAT.markers$cluster)){
cat('\n')
cat("### Cluster ",clust,"\n")
biomarkers <- getTopBiomarker(DAT.markers, clusterID=clust, topN=2)
plotBiomarkers(DAT, biomarkers, clust)
cat('\n')
} ##Construct the plot object
volcanoPlot <- function(DEG_df, caption="", topFC_labeled=5){
DEG_df$sig<- ifelse( DEG_df$p_val_adj<0.05 & DEG_df$avg_logFC<1.5, "p_val_adj<0.05",
ifelse( DEG_df$p_val_adj<0.05 & DEG_df$avg_logFC>1.5, "p_val_adj<0.05 & avg_logFC>1.5",
"p_val_adj>0.05"
))
DEG_df <- arrange(DEG_df, desc(sig))
yMax <- max(-log10(DEG_df$p_val_adj)) + max(-log10(DEG_df$p_val_adj))/3 #ifelse(max(-log10(DEG_df$p_val_adj))<45, 50, max(-log10(DEG_df$p_val_adj)) + 10)
vol <- ggplot(data=DEG_df, aes(x=avg_logFC, y= -log10(p_val_adj))) +
geom_point(alpha=0.5, size=3, aes(col=sig)) +
scale_color_manual(values=list("p_val_adj<0.05"="turquoise3",
"p_val_adj<0.05 & avg_logFC>1.5"="purple",
"p_val_adj>0.05" = "darkgray")) +
theme(legend.position = "none") +
xlab(expression(paste("Average ",log^{2},"(fold change)"))) +
ylab(expression(paste(-log^{10},"(p-value)"))) + xlim(-2,2) + ylim(0, yMax) +
## ggrepl labels
geom_text_repel(data= arrange(DEG_df, p_val_adj, desc(avg_logFC))[1:topFC_labeled,],
# filter(DEG_df, avg_logFC>=1.5)[1:10,],
aes(label=gene), color="black", alpha=.5,
segment.color="black", segment.alpha=.5
) +
# Lines
geom_vline(xintercept= -1.5,lty=4, lwd=.3, alpha=.5) +
geom_vline(xintercept= 1.5,lty=4, lwd=.3, alpha=.5) +
geom_hline(yintercept= -log10(0.05),lty=4, lwd=.3, alpha=.5) +
ggtitle(caption)
print(vol)
}
# Run plots
for (clust in unique(DAT.markers$cluster)){
cat('\n')
cat("### Cluster ",clust,": Volcano","\n")
cap <- paste("Cluster",clust,"DEG Table")
DEG_df <- subset(DAT.markers, cluster==as.character(clust)) %>% arrange(desc(avg_logFC))
volcanoPlot(DEG_df, caption = cap)
createDT(DEG_df, caption = cap)
cat('\n')
}##
## ### Cluster 0 : Volcano
##
##
## ### Cluster 1 : Volcano
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text_repel).
##
##
## ### Cluster 2 : Volcano
for (clust in top1$cluster){
subClust <- subset( top1, cluster==clust)
cat('\n')
cat("### Cluster",clust,"\n")
cat( "Biomarker\n",subClust$gene)
results <- Seurat::FindGeneTerms(QueryGene = subClust$gene)
print(results) #parse_html_notebook(results)
cat('\n')
}##
## ### Cluster 0
## Biomarker
## S100A12{xml_document}
## <html>
## [1] <body>\n<span class="term">CEBPA_26348894_ChIP-Seq_LIVER_Mouse</span ...
##
##
## ### Cluster 1
## Biomarker
## FCER1A{xml_document}
## <html>
## [1] <body>\n<span class="term">FOXA1_21572438_ChIP-Seq_LNCaP_Human</span ...
##
##
## ### Cluster 2
## Biomarker
## FCGR3A{xml_document}
## <html>
## [1] <body>\n<span class="term">TTF2_22483619_ChIP-Seq_HELA_Human</span> ...
fp <- FeaturePlot(object = DAT, features.plot = top1$gene, cols.use = c("grey", "purple"),
reduction.use = "tsne", nCol = nCols, do.return = T)top5 <- DAT.markers %>% group_by(cluster) %>% top_n(5, avg_logFC)
# setting slim.col.label to TRUE will print just the cluster IDS instead of
# every cell name
DoHeatmap(object = DAT, genes.use = top5$gene, slim.col.label=T, remove.key=T) # ifelse(interactive, p %>% ggplotly() %>% toWebGL() %>% print(), print(p))RidgePlot(DAT, features.plot = top1$gene, nCol = nCols, do.sort = F)## Picking joint bandwidth of 0.291
## Picking joint bandwidth of 0.13
## Picking joint bandwidth of 0.0842
Visualize biomarker expression for each cluster, by disease
top2 <- DAT.markers %>% group_by(cluster) %>% top_n(2, avg_logFC)
sdp <- SplitDotPlotGG(DAT, genes.plot = top2$gene, cols.use = c("blue","red"),
x.lab.rot = T, plot.legend = T, dot.scale = 8, do.return = T, grouping.var = "dx")The following plots show the absolute expression of each biomarker, as opposed to avg_logFC which is dependent on the expression patterns of other cell types being compared.
markerList <- c("CD14", "FCGR3A")
get_markerDF <- function(DAT, markerList, meta_vars =c("barcode", "dx", "mut","post_clustering", "percent.mito","nGene", "nUMI")){
exp <- DAT@scale.data %>% data.frame()
marker.matrix <- exp[row.names(exp) %in% markerList, ]
marker.matrix$Gene <- row.names(marker.matrix)
markerMelt <- reshape2:::melt.data.frame(marker.matrix, id.vars = "Gene", variable.name = "Cell",value.name = "Expression")
metaSelect <- DAT@meta.data[,meta_vars]
markerDF <- merge(markerMelt,metaSelect, by.x="Cell", by.y="barcode")
return(markerDF)
}
markerDF <- get_markerDF(DAT, markerList)
createDT(markerDF, caption = "Known Marker Expression")# Explore expression differences between groups
marker_vs_metadata <- function(markerDF, meta_var){
# Create title from ANOVA summary
ANOVAtitle <- function(markerDF, marker){
nTests <- length(unique(markerDF$Gene))
res <- anova(lm(data = subset(markerDF, Gene==marker),
formula = Expression ~ eval(parse(text=meta_var))))
title <-paste(paste("ANOVA (",marker, " vs. ",meta_var, ")", sep=""),
": p=",round(res$`Pr(>F)`,3),
", F=",round(res$`F value`,3),
ifelse(res$`Pr(>F)`<.05/nTests,"(Significant**)",
"(Non-significant)") )
}
title = ""
for (marker in unique(markerDF$Gene) ){
cat(marker)
title <- paste(title, "\n", ANOVAtitle(markerDF, marker))
}
ggplot(markerDF, aes(x=eval(parse(text=meta_var)), y=Expression, fill= Gene)) +
geom_violin() +
geom_point( position=position_jitterdodge(jitter.width = .2, dodge.width = .9 ), alpha=0.6, color="turquoise3") +
labs(title = title, x=meta_var) +
theme(plot.title = element_text( size=10)) +
scale_fill_manual(values=c("brown", "slategray"))
}marker_vs_metadata(markerDF, "dx")## CD14FCGR3A
marker_vs_metadata(markerDF, "mut")## CD14FCGR3A
identify_cellTypes_by_biomarkers <- function(DAT.markers, topN_search=5){
top <- DAT.markers %>% group_by(cluster) %>% top_n(topN_search, avg_logFC)
clust_cellTypes <- list()
for (clust in top$cluster){
clustSub <- top[top5$cluster==clust, ]
CD16_logFC <- subset(clustSub, gene=="CFD")$avg_logFC
cellType <- ifelse(sum(markerList %in% clustSub$gene), # Both CD14 and CD16? Great, keep going
ifelse(CD16_logFC == abs(CD16_logFC), "CD14++/CD16+", # But does CD16 have pos logfC? If so, then it's "CD14++/CD16+"
"CD14++/CD16--"), # Otherwise, it means it means CD16 logFC is neg, meaning "CD14++/CD16--"
NA) # If it's none of these, it's an undefined cell type
clust_cellTypes[clust] <- as.factor(cellType)
}
newMeta <- DAT@meta.data
newMeta["CellType_DGE"] <- plyr::mapvalues(newMeta$post_clustering, names(clust_cellTypes), as.character(clust_cellTypes) )
DAT <- AddMetaData(DAT, metadata = newMeta)
return(DAT)
}
DAT <- identify_cellTypes_by_biomarkers(DAT.markers, 5)
# (Doesn't make sense to do bar plot because whole clusters are defined by their biomarkers)
tSNE_metadata_plot("CellType_DGE")## t-SNE Metadata plot for CellType_DGE
# A simplistic way of categorizing cells into CD14++/CD16+ and CD14++/CD16--,
## is by splitting cells into groups based on whether their expression is
## higher or lower than the average CD16 expression of all cells.
identify_cellTypes_by_avgExpression <- function(DAT, markerDF){
avgMarkerExp <-markerDF %>% group_by(Gene) %>% dplyr::summarise(meanExp = mean(Expression))
avgMarkerExp <- setNames(avgMarkerExp$meanExp, avgMarkerExp$Gene)
CD16 <- markerDF[markerDF$Gene=="FCGR3A",]
CD16_group <- ifelse(CD16$Expression >= avgMarkerExp["FCGR3A"], "CD14++/CD16+", "CD14++/CD16--")
CD16["CellType_AvgExp"] <- as.factor(CD16_group)
# Make sure row order is same before putting back into meta.data
metaD <- DAT@meta.data
newMeta <- merge(metaD, CD16[,c("Cell","CellType_AvgExp")], by.x="barcode", by.y="Cell")
row.names(newMeta) <- row.names(metaD)
DAT <- AddMetaData(DAT, metadata = newMeta)
return(DAT)
}
DAT <- identify_cellTypes_by_avgExpression(DAT, markerDF)
# Get proportions of cell types in each cluster
cluster_proportions <- DAT@meta.data %>% group_by(CellType_AvgExp, post_clustering) %>%
tally() %>%
group_by(post_clustering, CellType_AvgExp) %>%
mutate(percentTotal = n/sum(n)*100)
ggplot(cluster_proportions, aes(x=post_clustering, y=percentTotal, fill=CellType_AvgExp)) + geom_col(position = "fill") +
labs(title="Proportions of Cell-types per Cluster: \n CellType_AvgExp",
x="Cluster", y="Cell Type / Total Cells") +
scale_fill_manual(values=c("brown", "slategray")) tSNE_metadata_plot("CellType_AvgExp")## t-SNE Metadata plot for CellType_AvgExp
markerDF <- get_markerDF(DAT, markerList,
meta_vars =c("barcode", "dx", "mut","ID","post_clustering", "percent.mito","nGene", "nUMI",
"CellType_DGE","CellType_AvgExp"))
Spectral <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(length(unique(DAT@meta.data$mut)), "Spectral"))
if (T==F){
# Spectral <- heatmaply::Spectral(length(unique(DAT@meta.data$mut)))
markerMelt <- reshape2::acast(markerDF, Cell~Gene, value.var="Expression", fun.aggregate = mean, drop = F, fill = 0)
heatmaply::heatmaply(markerMelt, key.title="Expression",#plot_method= "ggplot",
k_row = dim(markerMelt)[2], dendrogram = "row",
showticklabels = c(T, F), xlab = "Known Markers", ylab = "Cells", column_text_angle = 45,
row_side_colors = DAT@meta.data[,c("dx","mut", "CellType_DGE")], row_side_palette = Spectral
) %>% colorbar(tickfont = list(size = 12), titlefont = list(size = 14), which = 2) %>%
colorbar(tickfont = list(size = 12), titlefont = list(size = 14), which = 1)
}else{
# markerDF_sub <-subset(markerDF, Gene==markerList[1])
# var_to_colors(markerDF_sub, "post_clustering")
# library(pheatmap)
# pheatmap(markerMelt, annotation_row = markerDF_sub[c("dx","mut","CellType_DGE", "CellType_AvgExp")])
# pheatmap(markerMelt, kmeans_k = NA, annotation_row = markerDF_sub[c("dx","mut","CellType_DGE", "CellType_AvgExp")],
# cluster_cols = F, cutree_rows = length(unique(markerDF$post_clustering)), angle_col=45 )
library(RColorBrewer)
var_to_colors <- function(markerDF, metaVar){
colors <- brewer.pal(length(unique(markerDF[metaVar]) ), "Dark2")
sample(colors, length(unique(markerDF[metaVar])), replace = TRUE, prob = NULL)
# metaColors <- colors[ subset(markerDF, Gene==markerList[1])[metaVar][,1] %>% as.factor() ]
return(metaColors)
}
# library(GMD)
# myCols = cbind(var_to_colors(markerDF, "dx"), var_to_colors(markerDF, "mut"))
# rlab=t(cbind(
# var_to_colors(markerDF, "post_clustering"),
# var_to_colors(markerDF, "dx")
# ))
# heatmap.2(marker.matrix, key.title="Expression", col = viridis(300), trace="none",Colv = F, Rowv = F,
# labRow = F, xlab = "Biomarker", ylab="Cell", cexCol=1, RowSideColors = var_to_colors(markerDF, "post_clustering")
# )
# heatmap.3(markerMelt, dendrogram = 'row', kr = length(unique(markerDF)), labRow = F,
# xlab = "Biomarker", ylab = "Cell", RowSideColors = rlab, RowSideColorsSize=2 )
markerDF <- markerDF %>%
mutate_at(vars(post_clustering, dx, mut, ID, CellType_DGE, CellType_AvgExp), as.factor) %>%
mutate(Cluster = post_clustering) %>%
arrange(post_clustering)
# markerMelt <- reshape2::acast(markerDF, Cell~Gene, value.var="Expression", fun.aggregate = mean, drop = F, fill = 0)
markerMelt <- dcast(markerDF, Cell + post_clustering + dx + mut + ID + CellType_DGE + CellType_AvgExp ~ Gene,
fun.aggregate = mean, value.var = "Expression") %>% arrange(post_clustering)
marker.matrix <- markerMelt[markerList] %>%as.matrix()
row.names(marker.matrix) <- markerMelt$Cell
ha = HeatmapAnnotation(df = markerDF[c("dx","mut","ID","CellType_DGE","CellType_AvgExp","post_clustering")], which = "row")
ComplexHeatmap::Heatmap(marker.matrix, col=viridis(300), column_title = "Biomarker", row_title = "Cell",
row_dend_reorder = F,show_row_names = F, show_column_dend = F,show_row_dend =T,
cluster_rows = T, column_title_side = "bottom",km = length(unique(markerMelt$post_clustering))) + ha
} markerDF <- markerDF %>% mutate(Cluster = post_clustering)
# Show mean exp for each marker
avgMarker <- markerDF %>% group_by(Gene, Cluster) %>% summarise(meanExp = mean(Expression))
p <- ggplot(data = avgMarker, aes(x=Gene, y=Cluster, fill=meanExp)) %>% + geom_tile() + scale_fill_viridis()
p# ifelse(interactive, p %>% ggplotly() %>% toWebGL() %>% print(), print(p))# Show mean exp for each marker
avgMarker <- markerDF %>% group_by(Gene, dx, Cluster) %>% summarise(meanExp = mean(Expression))
p <- ggplot(data = avgMarker, aes(x=Gene, y=dx, fill=meanExp)) %>% + geom_tile() + scale_fill_viridis()
p# ifelse(interactive, p %>% ggplotly() %>% toWebGL() %>% print(), print(p))ggplot(data = markerDF, aes(x=Cluster, y=Expression, fill=Gene)) %>%
+ geom_boxplot(alpha=0.5) %>% + scale_fill_manual(values=c("purple", "turquoise")) #, results = 'hide', fig.show='hide'
expressionTSNE <- function(DAT, marker, colors=c("grey", "red")){
FeaturePlot(object = DAT, features.plot = marker, cols.use = colors,
reduction.use = "tsne", nCol=2, do.return = T, dark.theme = T)[[1]]
# p <- ifelse(interactive, p %>% ggplotly() %>% toWebGL(), print(p))
}
plot_grid(expressionTSNE(DAT, markerList[1]),
expressionTSNE(DAT, markerList[2], colors=c("grey", "green")))current.cluster.ids <- unique(DAT.markers$cluster) #c(0, 1, 2, 3, 4, 5, 6, 7)
top1 <- DAT.markers %>% group_by(cluster) %>% top_n(1, avg_logFC)
new.cluster.ids <- top1$gene #c("CD4 T cells", "CD14+ Monocytes", "B cells", "CD8 T cells", "FCGR3A+ Monocytes", "NK cells", "Dendritic cells", "Megakaryocytes")
DAT@ident <- plyr::mapvalues(x = DAT@ident, from = current.cluster.ids, to = new.cluster.ids)
TSNEPlot(object=DAT, do.label=T, pt.size=0.5, do.return=T) # ifelse(interactive, p %>% ggplotly() %>% toWebGL() %>% print(), print(p)) # Available DGE methods:
## "wilcox", "bimod", "roc", "t", "tobit", "poisson", "negbinom", "MAST", "DESeq2"
runDGE <- function(DAT, meta_var, group1, group2, test.use="wilcox"){
#print(paste("DGE_allCells",meta_var,sep="_"))
DAT <- SetAllIdent(DAT, id = meta_var)
DAT <- StashIdent(DAT, save.name = meta_var)
DEGs <- FindMarkers(DAT, ident.1=group1, ident.2=group2, test.use=test.use, only.pos = F)
DEGs$gene <- row.names(DEGs)
cap <- paste("DEGs:\n",group1, "vs.", group2)
createDT(DEG_df, caption = cap)
volcanoPlot(DEG_df, caption = cap)
DAT <- SetAllIdent(DAT, id = "post_clustering")
return(DEGs)
}DEG_df <-runDGE(DAT, "dx", group1 = "PD", group2="control")DEG_df <-runDGE(DAT, "mut", "LRRK2", "PD")## Warning: Removed 1 rows containing missing values (geom_hline).
DEG_df <-runDGE(DAT, "CellType_DGE", "CD14++/CD16+", "CD14++/CD16--") ## Error in WhichCells(object = object, ident = ident.1): Identity : CD14++/CD16+ not found.
DEG_df <-runDGE(DAT, "CellType_AvgExp", "CD14++/CD16+", "CD14++/CD16--") DGE_within_clusters <- function(DAT, meta_var, group1, group2){
for (clust in unique(DAT@meta.data$post_clustering)){
# Subset cells by cluster
DAT_clustSub <- Seurat::SubsetData(DAT, accept.value = clust, subset.raw = T)
cat('\n')
cat("### ",paste("Cluster ",clust,": ",group1," vs. ", group2, sep="") , "\n")
DEG_df <-runDGE(DAT_clustSub, meta_var, group1, group2 )
cat('\n')
}
}DGE_within_clusters(DAT, "dx", "PD", "control") ##
## ### Cluster 1: PD vs. control
## Warning: Removed 1 rows containing missing values (geom_hline).
##
##
## ### Cluster 0: PD vs. control
## Warning: Removed 1 rows containing missing values (geom_hline).
##
##
## ### Cluster 2: PD vs. control
## Warning: Removed 1 rows containing missing values (geom_hline).
DGE_within_clusters(DAT, "mut", "LRRK2", "PD")##
## ### Cluster 1: LRRK2 vs. PD
## Warning: Removed 1 rows containing missing values (geom_hline).
##
##
## ### Cluster 0: LRRK2 vs. PD
## Warning: Removed 1 rows containing missing values (geom_hline).
##
##
## ### Cluster 2: LRRK2 vs. PD
## Warning: Removed 1 rows containing missing values (geom_hline).
DGE_within_clusters(DAT, "CellType_DGE", "CD14++/CD16+", "CD14++/CD16--") ##
## ### Cluster 1: CD14++/CD16+ vs. CD14++/CD16--
## Error in WhichCells(object = object, ident = ident.1): Identity : CD14++/CD16+ not found.
DGE_within_clusters(DAT, "CellType_AvgExp", "CD14++/CD16+", "CD14++/CD16--")##
## ### Cluster 1: CD14++/CD16+ vs. CD14++/CD16--
## Warning: Removed 1 rows containing missing values (geom_hline).
##
##
## ### Cluster 0: CD14++/CD16+ vs. CD14++/CD16--
## Warning: Removed 1 rows containing missing values (geom_hline).
##
##
## ### Cluster 2: CD14++/CD16+ vs. CD14++/CD16--
## Warning: Removed 1 rows containing missing values (geom_hline).
If you perturb some of our parameter choices above (for example, setting resolution=0.8 or changing the number of PCs), you might see the CD4 T cells subdivide into two groups. You can explore this subdivision to find markers separating the two T cell subsets. However, before reclustering (which will overwrite object@ident), we can stash our renamed identities to be easily recovered later.
new_resolution <- 3.0
orig_resolution <- paste("resolution",params$resolution,sep="_")
DAT <- StashIdent(object = DAT, save.name = orig_resolution)
## Warning in BuildSNN(object = object, genes.use = genes.use, reduction.type
## = reduction.type, : Build parameters exactly match those of already
## computed and stored SNN. To force recalculation, set force.recalc to TRUE.
DAT <- FindClusters(object = DAT, reduction.type = "pca", dims.use = 1:10,
resolution = new_resolution, print.output = F)## 4 singletons identified. 33 final clusters.
DAT <- StashIdent(object = DAT, save.name = "resolution_3.0")
plot1 <- TSNEPlot(object = DAT, do.return = TRUE, no.legend = TRUE, do.label = TRUE, label.size=labSize)
plot2 <- TSNEPlot(object = DAT, do.return = TRUE, group.by = "post_clustering",
no.legend = TRUE, do.label = TRUE, label.size=labSize)
plot_grid(plot1, plot2)res3.0_markers <- FindAllMarkers(object = DAT, min.pct = 0.25, thresh.use = 0.25, only.pos = F, test.use = "wilcox")
top1_res3.0 <- res3.0_markers %>% group_by(cluster) %>% top_n(1, avg_logFC)
FeaturePlot(object = DAT, features.plot = top1_res3.0$gene, cols.use = c("green", "blue"))# Set back to orig
DAT <- SetAllIdent(object = DAT, id = orig_resolution) enrichr_dbs <- c("KEGG_2018", "Reactome_2016",
"GO_Biological_Process_2018", "GO_Molecular_Function_2018", "GO_Cellular_Component_2018",
"Rare_Diseases_AutoRIF_ARCHS4_Predictions", "Human_Gene_Atlas")
createDT(enrichR::listEnrichrDbs(), "Enrichr Databases")for (clust in unique(DAT.markers.sig$cluster)){
cat('\n')
cat("### Cluster ",clust,"{.tabset .tabset-fade}\n")
geneList <- subset(DAT.markers.sig, cluster==clust)$gene
results <- enrichr(genes = geneList, databases = enrichr_dbs )
for (db in enrichr_dbs){
cat('\n')
cat("#### ",db,"\n")
print(createDT(subset(results[[db]], Adjusted.P.value<=0.05), paste("Enrichr Results: ",db,"Cluster ", clust)))
cat('\n')
}
cat('\n')
}##
## ### Cluster 0 {.tabset .tabset-fade}
## Uploading data to Enrichr... Done.
## Querying KEGG_2018... Done.
## Querying Reactome_2016... Done.
## Querying GO_Biological_Process_2018... Done.
## Querying GO_Molecular_Function_2018... Done.
## Querying GO_Cellular_Component_2018... Done.
## Querying Rare_Diseases_AutoRIF_ARCHS4_Predictions... Done.
## Querying Human_Gene_Atlas... Done.
## Parsing results... Done.
##
## #### KEGG_2018
##
##
## #### Reactome_2016
##
##
## #### GO_Biological_Process_2018
##
##
## #### GO_Molecular_Function_2018
##
##
## #### GO_Cellular_Component_2018
##
##
## #### Rare_Diseases_AutoRIF_ARCHS4_Predictions
##
##
## #### Human_Gene_Atlas
##
##
##
## ### Cluster 1 {.tabset .tabset-fade}
## Uploading data to Enrichr... Done.
## Querying KEGG_2018... Done.
## Querying Reactome_2016... Done.
## Querying GO_Biological_Process_2018... Done.
## Querying GO_Molecular_Function_2018... Done.
## Querying GO_Cellular_Component_2018... Done.
## Querying Rare_Diseases_AutoRIF_ARCHS4_Predictions... Done.
## Querying Human_Gene_Atlas... Done.
## Parsing results... Done.
##
## #### KEGG_2018
##
##
## #### Reactome_2016
##
##
## #### GO_Biological_Process_2018
##
##
## #### GO_Molecular_Function_2018
##
##
## #### GO_Cellular_Component_2018
##
##
## #### Rare_Diseases_AutoRIF_ARCHS4_Predictions
##
##
## #### Human_Gene_Atlas
##
##
##
## ### Cluster 2 {.tabset .tabset-fade}
## Uploading data to Enrichr... Done.
## Querying KEGG_2018... Done.
## Querying Reactome_2016... Done.
## Querying GO_Biological_Process_2018... Done.
## Querying GO_Molecular_Function_2018... Done.
## Querying GO_Cellular_Component_2018... Done.
## Querying Rare_Diseases_AutoRIF_ARCHS4_Predictions... Done.
## Querying Human_Gene_Atlas... Done.
## Parsing results... Done.
##
## #### KEGG_2018
##
##
## #### Reactome_2016
##
##
## #### GO_Biological_Process_2018
##
##
## #### GO_Molecular_Function_2018
##
##
## #### GO_Cellular_Component_2018
##
##
## #### Rare_Diseases_AutoRIF_ARCHS4_Predictions
##
##
## #### Human_Gene_Atlas
eigengenes <- read.delim("Data/bulkMonocytes_WGCNAmodules_geneMembership.txt", row.names = NULL)
modules <- read.delim("Data/bulkMonocytes_WGCNAmodules_geneModules.txt", row.names = NULL, sep = "",
col.names = c("Ensembl","moduleColors"))
modules <- base::merge(eigengenes, modules,by="Ensembl" )
for (mod in unique(modules$moduleColors)){
cat('\n')
cat("### Module ",mod,"{.tabset .tabset-fade}\n")
geneList <- subset(modules, moduleColors==mod)$symbol
results <- enrichr(genes = geneList, databases = enrichr_dbs )
for (db in enrichr_dbs){
cat('\n')
cat("#### ",db,"\n")
print(createDT(subset(results[[db]], Adjusted.P.value<=0.05), paste("Enrichr Results:",db,"Module", mod)))
cat('\n')
}
cat('\n')
}##
## ### Module blue {.tabset .tabset-fade}
## Uploading data to Enrichr...
## Warning in enrichr(genes = geneList, databases = enrichr_dbs): genes must
## be a vector of gene names or a dataframe with genes and score.
## Done.
## Querying KEGG_2018... Done.
## Querying Reactome_2016... Done.
## Querying GO_Biological_Process_2018... Done.
## Querying GO_Molecular_Function_2018... Done.
## Querying GO_Cellular_Component_2018... Done.
## Querying Rare_Diseases_AutoRIF_ARCHS4_Predictions... Done.
## Querying Human_Gene_Atlas... Done.
## Parsing results... Done.
##
## #### KEGG_2018
##
##
## #### Reactome_2016
##
##
## #### GO_Biological_Process_2018
##
##
## #### GO_Molecular_Function_2018
##
##
## #### GO_Cellular_Component_2018
##
##
## #### Rare_Diseases_AutoRIF_ARCHS4_Predictions
##
##
## #### Human_Gene_Atlas
##
##
##
## ### Module red {.tabset .tabset-fade}
## Uploading data to Enrichr...
## Warning in enrichr(genes = geneList, databases = enrichr_dbs): genes must
## be a vector of gene names or a dataframe with genes and score.
## Done.
## Querying KEGG_2018... Done.
## Querying Reactome_2016... Done.
## Querying GO_Biological_Process_2018... Done.
## Querying GO_Molecular_Function_2018... Done.
## Querying GO_Cellular_Component_2018... Done.
## Querying Rare_Diseases_AutoRIF_ARCHS4_Predictions... Done.
## Querying Human_Gene_Atlas... Done.
## Parsing results... Done.
##
## #### KEGG_2018
##
##
## #### Reactome_2016
##
##
## #### GO_Biological_Process_2018
##
##
## #### GO_Molecular_Function_2018
##
##
## #### GO_Cellular_Component_2018
##
##
## #### Rare_Diseases_AutoRIF_ARCHS4_Predictions
##
##
## #### Human_Gene_Atlas
##
##
##
## ### Module turquoise {.tabset .tabset-fade}
## Uploading data to Enrichr...
## Warning in enrichr(genes = geneList, databases = enrichr_dbs): genes must
## be a vector of gene names or a dataframe with genes and score.
## Done.
## Querying KEGG_2018... Done.
## Querying Reactome_2016... Done.
## Querying GO_Biological_Process_2018... Done.
## Querying GO_Molecular_Function_2018... Done.
## Querying GO_Cellular_Component_2018... Done.
## Querying Rare_Diseases_AutoRIF_ARCHS4_Predictions... Done.
## Querying Human_Gene_Atlas... Done.
## Parsing results... Done.
##
## #### KEGG_2018
##
##
## #### Reactome_2016
##
##
## #### GO_Biological_Process_2018
##
##
## #### GO_Molecular_Function_2018
##
##
## #### GO_Cellular_Component_2018
##
##
## #### Rare_Diseases_AutoRIF_ARCHS4_Predictions
##
##
## #### Human_Gene_Atlas
##
##
##
## ### Module green {.tabset .tabset-fade}
## Uploading data to Enrichr...
## Warning in enrichr(genes = geneList, databases = enrichr_dbs): genes must
## be a vector of gene names or a dataframe with genes and score.
## Done.
## Querying KEGG_2018... Done.
## Querying Reactome_2016... Done.
## Querying GO_Biological_Process_2018... Done.
## Querying GO_Molecular_Function_2018... Done.
## Querying GO_Cellular_Component_2018... Done.
## Querying Rare_Diseases_AutoRIF_ARCHS4_Predictions... Done.
## Querying Human_Gene_Atlas... Done.
## Parsing results... Done.
##
## #### KEGG_2018
##
##
## #### Reactome_2016
##
##
## #### GO_Biological_Process_2018
##
##
## #### GO_Molecular_Function_2018
##
##
## #### GO_Cellular_Component_2018
##
##
## #### Rare_Diseases_AutoRIF_ARCHS4_Predictions
##
##
## #### Human_Gene_Atlas
##
##
##
## ### Module pink {.tabset .tabset-fade}
## Uploading data to Enrichr...
## Warning in enrichr(genes = geneList, databases = enrichr_dbs): genes must
## be a vector of gene names or a dataframe with genes and score.
## Done.
## Querying KEGG_2018... Done.
## Querying Reactome_2016... Done.
## Querying GO_Biological_Process_2018... Done.
## Querying GO_Molecular_Function_2018... Done.
## Querying GO_Cellular_Component_2018... Done.
## Querying Rare_Diseases_AutoRIF_ARCHS4_Predictions... Done.
## Querying Human_Gene_Atlas... Done.
## Parsing results... Done.
##
## #### KEGG_2018
##
##
## #### Reactome_2016
##
##
## #### GO_Biological_Process_2018
##
##
## #### GO_Molecular_Function_2018
##
##
## #### GO_Cellular_Component_2018
##
##
## #### Rare_Diseases_AutoRIF_ARCHS4_Predictions
##
##
## #### Human_Gene_Atlas
##
##
##
## ### Module black {.tabset .tabset-fade}
## Uploading data to Enrichr...
## Warning in enrichr(genes = geneList, databases = enrichr_dbs): genes must
## be a vector of gene names or a dataframe with genes and score.
## Done.
## Querying KEGG_2018... Done.
## Querying Reactome_2016... Done.
## Querying GO_Biological_Process_2018... Done.
## Querying GO_Molecular_Function_2018... Done.
## Querying GO_Cellular_Component_2018... Done.
## Querying Rare_Diseases_AutoRIF_ARCHS4_Predictions... Done.
## Querying Human_Gene_Atlas... Done.
## Parsing results... Done.
##
## #### KEGG_2018
##
##
## #### Reactome_2016
##
##
## #### GO_Biological_Process_2018
##
##
## #### GO_Molecular_Function_2018
##
##
## #### GO_Cellular_Component_2018
##
##
## #### Rare_Diseases_AutoRIF_ARCHS4_Predictions
##
##
## #### Human_Gene_Atlas
Determine whether each of the clusters in scRNA-seq data are enriched for WGCNA eigengenes (a weighted vector of all genes representing each co-expression module).
#Get the average expression of every gene in each cluster
allGenes <- get_markerDF(DAT, markerList = row.names(DAT@scale.data), meta_vars = c("post_clustering", "barcode") )
clusterGeneAvg <- allGenes %>% group_by(post_clustering, Gene) %>% summarise(meanExp = mean(Expression))
eigengenes_filt <- subset(eigengenes,symbol %in% unique(clusterGeneAvg$Gene))
clusts_by_mods <- base::merge(clusterGeneAvg[c("Gene","meanExp")], eigengenes_filt[c("symbol", modName)],
by.x="Gene", by.y="symbol")
cor.test()
corrplot()
heatmap.2
f <- function(module){
eigengene <- eigengenes[paste0("MM", mod)]
means <- tapply(eigengenes, DAT@meta.data$post_clustering, mean, na.rm = T)
return(means)
}
modules <- c("blue", "brown", "green", "turquoise", "yellow")
plotdat <- sapply(modules, f)
matplot(plotdat, col = modules, type = "l", lwd = 2, xaxt = "n", xlab = "Seurat Cluster",
ylab = "WGCNA Module Eigengene")
axis(1, at = 1:16, labels = 0:15)
matpoints(plotdat, col = modules, pch = 21)library(RRHO) #BiocManager::install("RRHO")
list.length <- 100
list.names <- paste('Gene',1:list.length, sep='')
gene.list1<- data.frame(list.names, sample(100))
gene.list2<- data.frame(list.names, sample(100))
for (clust in unique(DAT.markers.sig)){
# Compare each cluster
subClust <- subset(DAT.markers.sig, cluster==clust) %>% arrange(desc(avg_logFC))
for (mod in unique(modules$moduleColors)){
# Sort genes by module membership
modName <-paste("MM",mod,sep="")
subMod <- subset(modules, moduleColors==mod) %>% arrange(desc(eval(parse(text = modName))))
maxGenes <- min(length(subClust$gene), subMod$symbol) %>% as.numeric()
list1 <- subClust[1:maxGenes, c("gene","FC")] %>% dplyr::rename(value=FC)
list2 <- subMod[1:maxGenes, c("symbol",modName)] %>% dplyr::rename(gene=symbol, value=modName)
RRHO_path <-file.path("RRHO_results",paste(paste("Cluster",clust,sep=""),"vs",modName,sep="_"))
dir.create(RRHO_path,recursive = T, showWarnings = F)
RRHO_results <- RRHO(list1=list1, list2=list2,
labels = c(paste("Cluster",clust,sep="_"), paste("Module",mod,sep="_")),
plots = T, alternative = "enrichment", outputdir = RRHO_path, BY=TRUE
)
lattice::levelplot(RRHO_results$hypermat)
# Pval testing
pval.testing <- pvalRRHO(RRHO_results, 50)
pval.testing$pval
xs<- seq(0, 10, length=100)
plot(Vectorize(pval.testing$FUN.ecdf)(xs)~xs, xlab='-log(pvalue)', ylab='ECDF', type='S')
lattice::levelplot(RRHO_results$hypermat.by)
}
}# Save results for EACH run (in their respective subfolders)
saveRDS(DAT, file=file.path(params$resultsPath, "cd14-processed.rds") )